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Abstract. 

We explore the transition from order to chaos for the Bohmian trajectories of a 
simple quantum system corresponding to the superposition of three stationary states 
in a 2D harmonic well with incommensurable frequencies. We study in particular the 
role of nodal points in the transition to chaos. Our main findings are: a) A proof of the 
existence of bounded domains in configuration space which are devoid of nodal points, 
b) An analytical construction of formal series representing regular orbits in the central 
domain as well as a numerical investigation of its limits of applicability, c) A detailed 
exploration of the phase-space structure near the nodal point. In this exploration 
we use an adiabatic approximation and we draw the flow chart in a moving frame 
of reference centered at the nodal point. We demonstrate the existence of a saddle 
point (called X-point) in the vicinity of the nodal point which plays a key role in the 
manifestation of exponential sensitivity of the orbits. One of the invariant manifolds 
of the X-point continues as a spiral terminating at the nodal point. We find cases of 
Hopf bifurcation at the nodal point and explore the associated phase space structure of 
the nodal point - X-point complex. We finally demonstrate the mechanism by which 
this complex generates chaos. Numerical examples of this mechanism are given for 
particular chaotic orbits, and a comparison is made with previous related works in the 
literature. 
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1. Introduction 

The formulation of quantum mechanics based on Bohm's trajectories (de Broghe 1926, 
Bohm 1952a,b) has attracted considerable interest in recent years because it offers a 
powerful tool to visualize quantum processes in terms of quantum orbits. According to 
the Bohmian interpretation (see Bohm and Hiley 1993 and Holland 1993 for reviews) 
the particles are guided by the Schrodinger field via deterministic equations of motion. 
In this approach the moving particles, together with the Schrodinger field, form the 
basic ingredients of objective realilty at the quantum level. On the other hand, in 
orthodox quantum mechanics the orbits do not refer to deterministic particles' motions, 
but they represent the streamlines of the probability current j = {ilj*'Vip — ip'Vip*)/2i 
(we set h = 1). Thus they represent a 'Lagrangian' (Holland 2005) or 'hydrodynamical' 
(Madelung 1926) description of the quantum probabihty flow, while the Schrodinger 
equation yields the Eulerian description of the same flow. At any rate, the descriptive 
power of the Bohmian approach is independent of its ontological interpretation. In fact, 
applications of Bohm's trajectories in the literature have so far successfully addressed 
such basic quantum processes as the two-slit experiment (Phillipidis 1979), the spin 
measurement via Stcrn-Gcrlach devices (Dewdney et al. 1986), tunneling through 
potential barriers (Hirschfcldcr ct al. 1974, Skodje et al. 1989, Loprcorc and Wyatt 
1999), ballistic electron transport (Beenakker and van Houten 1991), superfluidity 
(Feynman Lectures, Feynman et al. 1963), etc. (see Wyatt 2005 for a review of 
applications of quantum dynamics with trajectories). 

In the present paper we focus on one particular aspect of Bohm's theory that 
refers to the distinction of the Bohmian trajectories into regular and chaotic. A number 
of studies in the hterature (e.g. Diirr et al. 1992, Faisal and Schwengelbeck 1995, 
Parmentcr and Valentine 1995, dc Polavieja 1996, Dewdney and Malik 1996, lacomclli 
and Pcttini 1996, Frisk 1997, Konkcl and Makowski 1998, Wu and Sprung 1999, 
Makowski et al. 2000, Gushing 2000, Falsaperla and Fontc 2003, de Sales and Florencio 
2003, Wisniacki and Pujals 2005, Valentini and Westman 2005, Efthymiopoulos and 
Contopoulos 2006) have so far converged to the conclusion that generic quantum 
systems of more than one degrees of freedom are characterized by the coexistence, in 
the configuration space, of both regular and chaotic orbits. Some estabhshed results 
regarding this distinction are: 

a) The regular or chaotic character of the quantum mechanical orbits does not 
necessarily correlate with the character of the classical orbits of the same system, i.e., 
there are examples of systems with classically regular and quantum mechanically chaotic 
orbits, or vice versa (see Efthymiopoulos and Contopoulos 2006 for a review). 

b) The emergence of chaos is associated with the existence of 'nodal points' in the 
configuration space, i.e., points of the configuration space at which the Schrodinger field 
becomes null. In some particular examples it has been possible to identify a mechanism 
by which the approach of an orbit near a nodal point introduces chaos (Makowski et al. 
2000, Wisniaski and Pujals 2005). However the general problem of the mechanism of 
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transition from order to chaos for the Bohmian orbits is still largely unexplored. 

c) In some cases it has been shown that the extent of chaos is related to the 
number and spatial distribution of the nodal points in the configuration space (Frisk 
1997, Wisniacki et al. 2006). However, in other cases we find the manifestation of strong 
chaos even if only one nodal point is present. This problem is important because it has 
been shown that, when the degree of chaos is large, it is possible to obtain an asymptotic 
convergence of the distribution p of an ensemble of Bohmian trajectories to Born's rule 
p— I'^p, via a Bohm-Vigier (1954) stochastic mechanism (Valentini and Westman 2005, 
Efthymiopoulos and Contopoulos 2006), without having to postulate this rule. 

Our purpose in the present paper is to study how the transition from regular 
to chaotic motion takes place by considering the orbits of a very simple quantum 
system, namely the superposition of three eigenstates in a Hamiltonian system of two 
independent harmonic oscillators. Parmenter and Valentine (1995,1996) demonstrated 
that when the two oscillator frequencies are incommensurable, the configuration space 
is filled by both regular and chaotic orbits (the initial claim that all the orbits were 
chaotic (Parmenter and Valentine 1995) was corrected by Parmenter and Valentine 
(1996), and by Efthymiopoulos and Contopoulos (2006)). The same coexistence was 
found by Wisniacki and Pujals (2005) when the frequencies are commensurable but the 
ampitudes of the superposed eigenfunctions have a complex ratio, and by Konkel and 
Makowski (1998) in the case of a particle in a box with infinite walls. In the above cases 
there is only one nodal point which infiuenccs the orbits and introduces chaos. Our aim 
below is, then, to study the transition from order to chaos from two different points of 
view: a) topologically, i.e. we seek to distinguigh which domains of initial conditions 
lead to regular or chaotic motion and what theoretical criteria can be devised in order 
to separate these domains, and b) dynamically, i.e., we seek to identify the dynamical 
mechanism behind the transition from order to chaos. 

The following is an outline of the paper and of the main results: 

a) Section 2 contains a proof of the existence of domains in configuration space 
where nodal points cannot appear. In our considered example the size of these domains 
depends on the relative amplitudes of the three eigenfunctions and specific quantitative 
estimates of this dependence are given, which are compared to numerical results. 

b) In their domain of analyticity (i.e. far from nodal points) , the equations of motion 
admit solutions expandable in series of a properly defined small parameter (section 3). 
The series' terms can be determined by an iterative algorithm. The resulting solutions 
define theoretical orbits which are, by definition, regular. The theoretical orbits explain 
all the basic characteristics of the regular orbits as found by numerical integration. In 
particular they explain the frequencies, the form, the limits and the inner defiections of 
regular orbits. 

c) Section 4 passes to the other limit, of motion, close to the nodal points. In order 
to unravel the mechanism by which the orbits approaching the nodal point become 
chaotic, the key point is to take into account the motion of the nodal point itself by 
passing to a description of the orbits in a moving frame of reference centered at the 
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(moving) nodal point. The main characteristics of the orbits in this frame are found 
by expansions of the equations of motion in terms of a new small parameter, i.e., the 
distance R = e from the nodal point. The angular frequency of motion near the nodal 
point is of order 0(l/e^), a fact allowing us to use an adiabatic approximation. In this 
approximation, the flow lines near the nodal point are spirals terminating at the nodal 
point. However, the flow further away from the nodal point is quite complicated and 
it is studied in detail. In particular, we demonstrate how this flow is related to the 
manifestation of chaos in the system. The main results are derived theoretically and 
then substantiated by detailed numerical experiments. 

d) Finally, we discuss (section 5) how do our results compare with previous works 
in the literature on nodal points, and we end with the conclusions (section 6). 



2. Limits of nodal lines 

We study the quantum orbits in the Hamiltonian model of two uncoupled oscillators: 

H=l{pl+pl) + l{x' + {cyf) (1) 

when the guiding field is the superposition of three stationary states (Parmenter and 
Valentine 1995): 

V'(x, y, t) = e" A^-^^(l + axe-'' + bc'/'^xye-'^^^"^') . (2) 
The equations of motion (de Broglie 1926, Bohm 1952a) are 

(i,V}^^(^) (3) 



where 9 denotes the imaginary part of each of the components of the vector Vtp/tp, or 
dx a sint + bc^^^y sm{l + c)t 



dt G 
dy bc^^'^x{ax sin ct + sin(l + c)t) 

G 



(4) 



with 



G = l+2ax cost+2bc^^'^xy cos{l+c)t+a^x^+2abc^^'^x^y cosct+b'^cx^y^ .(5) 



The equations of motion (jl]) become singular whenever G = 0. From we then find 
the equations of the nodal points: 

sin(l + c)t asint 
asinct ' 6c^/^ sin(l + c)t 

When the frequency c is a rational number, the nodal points describe periodic motions in 
the configuration space {x, y), along a finite number of nodal lines (Figure la). However, 
when c is irrational there is an infinite number of nodal lines that fill open domains of 
the space {x,y) (Figure lb). As shown in section 3, a theoretical approximation of the 
regular orbits can be obtained in the complement of the domain of nodal lines. We thus 




Figure 1. Nodal lines given by Eqs.® for a = 1, = 1, < i < 1000 and (a) 
c = 7/10, (b) c = y2/2. 



first provide, in this section, rigorous bounds for tlie domain of nodal lines. In particular, 
we prove the following proposition: when the relative amplitudes a, h are non-zero, the 
domain of nodal lines is hounded by a set of limiting hyperbolae. 
To this end, we write Xq = —Xq/u and yo = —aYQ/bc^^"^ , where 

_ sin(l + c)t _ sint 

smcr sin(l + cji 

We make use of the following lemma: 

2 

(i) V n G 7^ with < m < 7r/2, 3 6 >0 with - < 6 < 1 : sin u = 6u 

TC 

(ii) V u eTZ with < m < tt, sinw <u. (8) 

In order to find the bounds of (|Xo|, |lo|) we consider all the 4 x 4 = 16 possible 
combinations of the trigonometric arguments t, {l + c)t being in any of the four quartiles 
of the trigonometric circle, i.e., 

(1 + c)t = 2kiTC ± ti, or (1 + c)t = {2ki + 1)% ± ti 

t = 2k2-n ± ta, or t = {2k2 + l)n ± (9) 

with [ki, k2) G Z"^ and < ti < 7r/2, < ta ^ 7i"/2. Taking now ct = (1 + c)t — t we 
find the following possibilities for | sin(ct)|: 

|sin(ct)| = |sin(ti ±^2)1 • (10) 

We then distinguish four cases, namely cases A,B with | sinct| = | sin(ti — ^2)!; and cases 
C,D with I sinct| = | sin(ti +^2)1- 

Case A: ti > t2 (and | sin(ct)| = sin(ti — 12)). Since ti — ^2 < 7i"/2, according to the 
Lemma ©(i) 

2 

3 (5i > 0, i = 1, 2, 3 : - < 5^ < 1 such that 

TT 
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(a) (b) 




Figure 2. The limits of nodal lines as determined by (a) Eq. p^ . (b) Eq. fT^ . (c) 
Eq.dlll), (d) Eqs.lHSl) and (e) The nodal lines given by Eq.® for a = 1, 6 = 1, 

c = V^/S, < i < 1000, superposed to the less restrictive of all the limits shown in 
panels (a) to (d). 
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I sin(l + c)t\ = 6iti, I sint| = ^2^2, I sinct| = sin(ti — ^2) = ^siti ~ ^2) 



or 



From (fTT!) then follows the inequality: 

which is shown graphically in Fig. 2a. When |lo| = 0, |Xo| is between 2/7r and 7t/2. 

Case B: ^2 > (and | sin(ct)| = sin(t2 — ^i))- Working in the same way as for Case 
A we find the inequality 

^ < |Xo| < 2,,,! , (13) 



fl>^o|-l-' ^'-||>o|-l 
shown graphically in Fig. 2b. When |Xo| 00, |yo| is between 2/7r and 7r/2 (dashed 
asymptotic curves in Fig. 2b). 

Case C: ti + ^2 < 7r/2 (and | sin(ct)| = sin(ti + ^2))- In this case we find that 

2 

3 5. > 0, i = 1, 2, 3 : - <5^<l such that 

vr 

I sin(l + c)t\ = 6iti, I sint| = 52^2, | sinct| = sin(ti + ^2) = 53{ti + t2) 



or 



^ < |Xo| < TTTRT^^ (14) 



TTdFol + l)-' ^'-2(|ro| + l) 
(Fig.2c). When |Fo| = 0, |Xo| is between 2/7r and 7r/2. 

Case D: ti + t2 > vr/2 (and | sin(ct)| = sin(ti + ^2))- In this case 

2 

3 5. > 0, z = 1, 2, 3 : - <5^<l such that 

TT 

I sin(H-c)t| = 6iti, I sint| = ^2^25 I sinct| = sin(ti+t2) = sin(7r — (ti+t2)) = 53(7r— ti— 12) 
implying 

2 1 
7r(2|Vo| + l) ^ - l^ol 

Furthermore, if ti < 7r/4 (case Di) we have ^2 > 7r/4 and vr — (ti + t2) > 7r/4. The 
Eqs.(I7]) and ([9]) imply that |Xo| < 1 and |Fo| > 1- On the other hand, if ti > 7r/4 (case 
D2) we have 

sint 1 

The union of the two possibilities yields: 

\Xo\ < 1 or iFol < 2 . (16) 

The inequalities (fT5|) and (fT6l) are shown graphically in Fig. (2d). When \Yq\ = 0, |Xo| 
is between 2/7r and 7r/2. 
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The less restrictive of all bounds determine the permissible domain for nodal points. 
The outer limit (|T3l) of Case B is less restrictive than the limit |Xo| of Case D if 
l/(^|^o| ~ 1) < 1; i-G- if l^ol < (in the exceptional case that IIqI > the limit 
|Xo| = 1 is less restrictive than the limit of case B). Figure 2e shows the analytical 
estimates for the bounds of the nodal lines compared to the numerical determination of 
the nodal lines in the variables (xo,i/o) for the parameters a = 6 = 1, c = ^2/2. The 
further restrictions of the cases A-D are also satisfied if ti and ^2 take the values specified 
in these particular cases. The main conclusion is that there is a central domain, with 
boundary defined by the innermost arcs of limiting hyperbolae, that is never crossed by 
nodal points. There is also an outer domain, at large distances from the center, which 
is again prohibited to nodal points. By calculating many orbits, our numerical evidence 
is that the orbits lying within these domains are regular. In particular, we now turn our 
attention to the analytical determination of the regular orbits for the central domain 
free of nodal points. 

3. Integrals of motion and the bounds of regular orbits 

According to the analysis of the previous section, a lower bound for the distance of 
a nodal point (xo,i/o) from the origin (0,0) is given by d^in = (a^o.mm + ?/o,mm)^''^ 
with (xo = Xo^rnin/a,yo = aYo,min/ hc^^'^) , and Xo,„i„,Fo,mm corresponding to the closest 
approach of the innermost limiting hyperbola (Eq.f llSI) ) to the origin, i.e. 

'^^-^0,min , T^XQ^min t _ n, \r _ ( ^ ^\ 

Ta ^ o - U, o,mm - - — i 

ib 2 TT \7rAo,mm / 

or Xo^mm = 0.461226, Yo^rnin = 0.242092. When a and b are large, the distance dmin 
is small and the inner domain free of nodal points is small. On the other hand, in the 
limit a — > 0, 6 we have dmin oo and the whole space is free of nodal points. 
This implies that if we look for analytic solutions of the equations of motion (jl]), i.e., 
free of singularities in the neighborhood of the origin, we may consider a,b in Eq.(l5]) 
as small parameters, and expand l/G in Eq.(jl]) in a power series of these parameters. 
This results in solutions 

x{t) =Xo + Xi{t)+X2{t) + y{t) = yo + yi{t)+y2{t) + ..., (17) 

where the functions x„(t), yn{t) are of order n in the amplitudes a or b. The convergence 
of these series, which are of the form of the 'third integral' (Contopoulos 1960), is an open 
problem. However, our numerical evidence below is that the form of theoretical orbits 
derived by these series fits well the form of the numerical orbits for small parameters 
a, b. 

In the zeroth order approximation all the solutions are equilibria x{t) = xq, 
y(t) = yo- This corresponds to the limit a = 6 = 0, in which the guiding t/^— field 
is a bound stationary state, and all the quantum orbits are neutral equilibria. On the 
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other hand, when higher order terms are taken into account, Eqs. fll7p can be inverted, 
and, provided that the series converge, they yield the integrals of motion 

X - xi{t) - X2{t) - ... = xo, y - yi{t) - y2{t) - ... = yo (18) 

implying that the resulting orbits are, by definition, regular. 

The solutions (fT7|) are found recursively, i.e., order by order, giving Xk, yu as exphcit 
trigonometric expressions in t, with two basic frequencies Ui = 1 and 1^2 = c. 

The first order equations read 

^ = -asint-6c^/2^osin((l + c)t), ^ = -60^/^X0 sin((l + c)t) (19) 
dt dt 

and they can be readily integrated yielding 

, , , , 1 /n cos((l + c)t) 
xi(t) = Ai + acost + &c^/ 1/0 — 

+ (20) 

where Ai,Bi are integration constants. Since we wish to absorb the whole dependence 
of the solution on the initial conditions in the Xo,yo part of the solution, we select the 
values of Ai, Bi such that Xi(0) = yi{0) = 0. 

In a similar way we treat the second order equations, finding the solutions 

CL^X b'^cx 
X2{t) = _ ^ cos(2t) - cos((l + c)t) 

2 (1 + c)^ 

b'^cxo / 2 1 

" 2(1 + c) " 20 + c 

2abc^^'^xoyo 



^°"2(TT^J'°'^^^^ + ''^^^ ^^^^ 

cos((2 + c)t) + A2 



and 



2 + c 

( abc^^"^ b'^cyo 



, . f abc^''' b^'cyQ \ 



X , a^n - r cos(2(l + c)t) (22) 

2(1 + c) V ° 2(1 + c)/ ^ ^ ' ' ^ ' 

ab , , abc^^"^ / n 1\ , , „ 

+ —j^ cos(ct) - —— - - cos((2 + c)t) + B2 



2ci/2 ^ ^ 2 + c 

respectively. The integration constants ^2,-82 are also given values such that ^2(0) = 

2/2(0) = 0. 

We can prove the consistency of the above construction, i.e., that no secular terms 
can appear in the above recursive scheme. The proof follows by induction: If the 
solutions Xi{t),yi{t), i = I, n contain only cosine terms (of the form cos((mi +m2c)t) 
with nil, 1712 integer, then, the expansion 



^ = i + E(-i)'= 
k=i 



2aa;(") cost + 26c^/23,(n)^(n) ^^^^^ ^ 
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(a) (b) (c) 




XXX 

Figure 3. Three orbits in the equations of motion ^ ioi a — h = 1, c = \/2/2 
and < < < 1000. The initial conditions are (a) a;(0) = 0.75, ?/(0) = 0.25 (regular, 
not overlapping with the domain of nodal lines), (b) x{Q) — y{0) — 1 (regular, partly 
overlapping with the domain of nodal lines), and (c) a;(0) = y(0) = 1.4 (chaotic). 



1 

+a2(x("))2 + 2a6c^/2^x("))2y(") cosct + b^c{x^''^)\y^''^)^ +... 

with 

= Mt) + xi{t) + ... + Xn{t), y^''\t) = yo{t) + yi{t) + ...+ y^{t) 

contains only powers of cosine terms, yielding again only cosine terms since 
COSW1COSW2 = (cos(wi + W2) + cos(wi — W2))/2. Thus, the equations of motion in 
order n + 1 yield only sine terms, since 1/G in Eqs.dlj) is multiplied only by sine 
terms and sinwi cos 1172 = (sin(wi + W2) + sin(wi — W2))/2. If for some terms we have 
Wi = (mi + m2c)t, W2 = {m[ + m2c)t, the sine terms produced in the equations of 
motion are sin(mi + m[ + m2C + mgc)^, or sin(mi — m[ + m2C — mgc)^. If (mi = ±m[ 
and m2 = ±m'2), or, c is rational and equal to c = —(mi ± m[)/{m2 ± iTL'2), one of the 
sine terms in the equation of motion becomes equal to zero (resonance). However, the 
resonances do not create any secular term in the solutions of the equations since the 
associated terms simply disappear from the r.h.s of Eqs.(jl]). Thus, the equations in the 
next order give dxn+i/dt and dyn+i/dt as sums of sine terms, and it follows that, if x^"^ 
and y^'^^ are sums of cosine terms, Xn+i and yn+i are also sums of cosine terms. Thus, 
no secular terms appear in the solutions x*-"-*, y^'^\ Vn = 1, . . . 00 since x^^\y^^^ are sums 
of cosine terms. 

The main characteristics of the regular orbits in the central region can be 
understood in terms of the above equations. In particular: 

a) The regular orbits are quasi-periodic, i.e., they are given as double Fourier series 
with two fundamental frequencies z/i, z/2, which have constant values i^i = 1, z/2 = c. 
This fact is important because it implies that there is no dependence of the frequencies 
of the orbits on the amplitudes of the oscillations. This is different from what is usually 
encountered in the case of classical nonlinear dynamical systems, in which the frequencies 
depend, in general, on the amplitudes. On the contrary, the quantum mechanical orbits 




I, iiyn I I — II I i J 1 I 1 I ' I (I'll! ^ ^ -Vi i^n f~r^ I J I I I II ' ' I I I I I ' ^ ^^1 I -Wi I l i i ^ i i i i i i i . ■ i i i i l n i i . . i i i m i i i i i ; 

X X t 

Figure 4. (a) An orbit with a = b = 0.2, c = \/2/2 and initial conditions a:(0) = 1, 
y{0) = 0. (b) The theoretical approximation to the same orbit up to terms of second 
degree in a and &, given by Eqs. ((23)) and (|24| . (c) The time evolution of y{t) for the 
same orbit. The numerical and theoretical curves almost coincide, (d) Same as (a) 
but with initial conditions a;(0) ~ y{0) — 0. (e) Same as (b) but for the orbit (d). (f) 
Same as (c) but for the orbit (d); the secondary local maxima of y{t) correspond to 
deflections of the orbit inside its 'box' limit. 



in this system are completely degenerate with respect to their fundamental frequencies. 

b) The amplitudes of all the trigonometric terms depend on a and b, but in many 
terms they depend also on the initial conditions XQ,yQ. The former condition ensures 
that for a, b sufficiently small a regular orbit does not overlap with the domain of nodal 
lines. This is because the amplitudes of oscillations are 0(a) and 0(6c^/^) in the x and 
y axes respectively (Eqs. (l20l) ). while the minimum distance of a nodal point from the 
center is of order 0(|l/a| + \l/b\) (section 2). Numerically we find such regular orbits 
for a,b as high as a = 6 = 1 (Figure 3a). However, the oposite is not true, namely an 
orbit overlaping partly with the domain of nodal lines may still be regular (Figure 3b). 
In fact, we find numerically that while there is a spatial overlap of the domains of the 
orbit and of the nodal lines, the time evolution of both the orbit and the nodal point 
is such that their distance is always large (of order unity). Thus when an orbit enters 
some region containing nodal lines the nodal point is far from this region. However, we 
also find numerically that if an orbit has significant overlap with the domain of nodal 
lines, the orbit is, in general, chaotic (Figure 3c). 



Nodal points and the transition from ordered to chaotic Bohmian trajectories 



12 



0.B9 




-0.30 



-0.50 



0.50 



(b) 



0.30- 



-0.30- 



.50- 




-1 .20 



-0.80 



-0.40 



0.00 



- 1 .20 



III I I I I r-r 

-0.80 -0.40 0.30 




0.60 



0.50 




-0.S0 



Figure 5. Numerical versus tlieoretical orbits for a = b = 0.5, c = \/2/2. 
(a) Numerical orbit with initial conditions x(0) = y(0) — 0. (b) The theoretical 
approximation of orbit (a) with an expansion up to terms of 10-th degree in a,&. (c) 
Numerical orbit with initial conditions a;(0) = 0.5, ?/(0) = 0. (d) Same as in (b) but 
for the orbit (c). 



c) In the lowest approximation, when 7^ the theoretical orbits are 'box orbits' 
(Eqs.( l20l) ) like in the classical case. However, when higher order terms are taken into 
account, some box orbits develop deflections internal to the box, while other box orbits 
have deflections only at the boxes' limits, depending on the value of xo,|/o- A simple 
example is provided by orbits starting on the axis Uq = 0. Then, the equations of the 
orbits are simplifled considerably . We have, up to second order: 

x{t) = Xo + a[cost — 1] — [cos2i(: — 1] — 
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''^^ :[cos((l + c)t) - 1] + --^^[cos(2(l + c)t) - 1] (23) 



[l + cy vv ■ . / J ■ 4(1 + c) 



^W = 3^[eos((l + i)t)-l] 



^ ^ ^ ,cos((l + c)t) - 1] + ^[cos(ct) - 1] (24) 

-^(^g-^)[cos((2 + c)t)-l] . 

In Eg. (123 P the first order term a [cost — 1] depends only on a, yielding an oscillation of 
amplitude 2a, while in Eq.( l24l) . the first order term 6c^''^xo[cos((l + c)t) — 1]/(1 + c) 
depends on both h and Xq. For all initial conditions with |xo| = 0(1), this term 
is dominant over the second order terms in Eq. fl2^ . Thus the orbit is a deformed 
parallelogram, i.e., it resembles to classical box orbits (Figure 4a,b,c). On the other 
hand, if |xo| —^0, this term becomes small, while the second order terms, not depending 
on xo, become now important. This means that secondary oscillations of the function 
y{t) are developed, that may also yield local minima or maxima besides the main 
minimum or maximum defined by the first order Fourier component of flM|) . This 
causes the orbit to develop deflections in the y— direction inside the box. A numerical 
example of this behavior is demonstrated in Figs. 4d,e,f. 

Figure 5 shows some further examples of theoretical orbits calculated by the above 
series, via a computer program implementing the recursive algorithm up to the 10th 
order. As expected in any kind of perturbative series, as a, h increase one needs higher 
order terms to obtain a good approximation of the orbits. In fact, as already analyzed, 
the approximation depends also on the values of the initial conditions Xq, Vq which appear 
in the amplitudes of the trigonometric terms of Eqs. fl20l) . (I2T1) and (1221) . Thus, when 
^0 = Ho = 0, the theoretical orbit for a = b = 0.5 (Fig.5a,b) fits well the numerical 
orbit, while if the series are truncated at orders well below 10 the agreement is not 
so good. On the other hand for the same amplitudes a, b the fit is not good when 
2^0 = Z/o = 0.5 (Fig. 5c, d). Now, in any perturbation theory the analytical results are 
precise up to values of the small parameters smaller than the values for which regular 
orbits are found numerically. Thus, in our case we find that for amplitudes larger than 
a = 6 = 0.75 the approximation of the numerical orbits by series yields no longer 
accurate results, despite the fact that we still find numerically many regular orbits. 



4. The dynamics close to nodal points and the transition to chaos 

4.I. Phase space structure close to nodal points and the adiabatic approximation 

Having demonstrated the existence of regular orbits far from the nodal points, we now 
examine the motion in the other limit, i.e., close to a nodal point. Our main remark 
in the sequel is that the phenomena relevant to the transition to chaos are unraveled 
when one considers the passage of the orbits in the neighborhood of the nodal point in 



Nodal points and the transition from ordered to chaotic Bohniian trajectories 



14 



a moving frame of reference that is centered at the nodal point. Introducing 
V = y — yo, the equations of motion in the moving frame read: 

du hc^/'^v sin(l + c)t 

~dt^ G 

dv hc^^'^u sin(l + c)t — abc'^/'^v? sin ct 

Tt = G ^25) 

where G = G2 + G^ + G^ with 

= — - 2bc^^'^uv cos(l + c)t + fe^cXgW^ 

Xq 

G3 = u^v cos(l + c)t + 2h'^cxQUv'^ (26) 

and xo,2/o are found by differentiating xq, yo with respect to time in (Eq.(l6]). We then 
consider the distance R = \Ju^ ^ = e of an orbit from the nodal point as a small 
parameter and derive the main characteristics of the motion by taking expansions of the 
equations of motion (|25l) in powers of e, i.e., by considering both u and v as 0{K) = 0(e). 
In polar coordinates u = R cos (j), v = R sin (p, the equations fl25|) read: 

dR ahc^l^R^ cos^ sin sin ct 

— = Xo cos - ?/o sm 

at G 

(i0 foc-*^/^ sin(l + c)t — afec-'^/^-Rcos^ 0sinct 1 1 

rft = Q -^y^ cos ^ + ;^^o sm (27) 

where G = g2R^ + gsR^ + g^iR'^ and g2,g3,g4 are readily specified from Eq.fl26l) (see 
appendix for explicit formulae). 

To the leading order (1/e^), the second of Eqs. (!27l) yields: 

# bc'/^ sin(l + c)t 

17 ~ 7 \" ■■■ ^ ) 

dt ^2 1 _ 25ci/2 cos sin cos(l + c)t + h^cxl sin^ j 

Thus, the angular velocity around the nodal point is large (of order 0(l/e^)) and the 
period can be made arbitrarily small by approaching closer and closer to the nodal 
point. This fact justifies the use of the adiabatic approximation in the study of the 
motions near the nodal point. That is, at a given initial time to we set t = to + e^t', 
with e = Ro = R(to) , and find 

d(j) bc^ sm{l + c)to 



df ^ - 26cV2 cos sin cos(l + c)to + b'^cxl sin^ 



+ 0{e') . (29) 

Now, the denominator in the r.h.s. of Eq. fl29|) is equal to the square of the length of 
one diagonal of the parallelepiped with sides | cos0/xo|,|&c^''^xo sin 0| forming an angle 
(1 + c)to, thus it is always positive. It follows that d(j)/dt' and dcp/dt have a unique sign 
during a whole period of 0, which is the same as the sign of sin((l + c)to. That is, at a 
given time to-, the angular motions close to the nodal point are all described in the same 
sense. 
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In the same approximation we can now determine the form of the integral curves of 
the velocity vector field fl25l) . Dividing the first with the second of Eqs. fl27|) . and setting 
a constant to in the place of t in the r.h.s. yields the equation of the integral curves, 
which is of the form: 

dR ^ ^2(0; to, xq, xq, yo)R^ + ^3(0; tp, xp, xp, y^jR^ + ■■■ , 
dcj) 5p(0;tp,xo,xo,?/o) + 5i(0;tp,Xp,Xp,?/p)/?+ ... 
The precise functions Aj, Bi are given in the appendix. At this point it suffices to state 
that both functions contain only trigonometric terms in (p. If we expand Eq. (|30l) with 
respect to R we find: 
dR 

-J- = /2(0; to, a;p, xp, yQ)R^ + h{(t)] tp, xp, xp, ?/p)i?^ + Oi^R^) (31) 
d(p 

with 

i?o' 5o Bl ■ 

Rescaling the radial distance as i? = ei?', with e = i?p, Eq. fl3Tl) takes the form: 



e/2/?'' + eV3i?" + 0(e^i?'^) . (32) 

This equation satisfies the neccesary conditions for applying the averaging theorem (e.g. 
Verhulst 1993). This implies that there is a near-identity transformation R' = R' + 0{e) 
such that the dynamics in terms of R' is given by: 

dP' 

^ = e<f2>R" + e^<h>R'^ + Oie^R'^) (33) 
a(p 

where 



1 

< fi>=7r fi^^^ i = 2,3 

ZTT JO 



After some algebra (see appendix) we find that < /2 >= and < /a >^ 0. Thus, 
the equation of the integral curves (back-transformed to non-rescaled variables) reads 
finally: 

^=<h>R' + ... (34) 

where < /s > depends only on tp both exphcitly and through xp,xp,?/p. The precise 
form of < /s >, found in the appendix, reads: 

^ / 1 + b^cx^o \ 
l^46ci/2x4sin(l + c)tpyl ^ 

/ . , a;p?/p(6^cxQ — 1) 2/ -2 ■2\ , \ . i /qp-n 

r "° + 6cV2sin(l + c)tp - - + ) • ^''^ 



Eq. (lMl) admits the solution 

i?(0) = , ^° (36) 

y'l-2i?2</3>(0_0o) 

which is a spiral terminating at i? = 0, i.e., at the nodal point, when — >■ oo (if 
< fs > < 0), or ^ — oo (if < fs > > 0). Comparing this to the actual sense of 
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rotation, given by the sign of sin((l + c)to, we can decide whether, as t increases, the 
orbits along the spiral recede from or approach to the nodal point. The sense of rotation 
changes at times to when sin(l + c)to = 0, i.e., xq = and = oo. On the other hand, 
the sign of < /s > changes at times to when 

^oxo + f T^^, u ~ ^o(io - Vo) cot(l + c)to = . 
bc^'^ sm(l + cjto 

Close to such a time a Hopf bifurcation takes place that is connected to a change of the 
character of the nodal point from attractor to repellor, or vice versa. This phenomenon 
is analyzed in subsection 4.2 below. 

Eg. (1361) is valid only very close to the nodal point. In order to find the form of the 
phase flow at larger distances from the nodal point, we look for stationary points of the 
flow (1251) . The stationary points are given by non-zero solutions {uq, Vq) of the system of 
equations du/dt = dv/dt = 0. Assuming {uo,Vo) small, and keeping terms up to second 
degree in Uq, Vq in Eqs. fl25l) . we find the solution: 

[ io \ [ ct'^o '^^0 \ 
^° ~ \yoJ \sm{l + c)to ^7 

which, after replacement in the first of Eqs. (1251) . with du/dt = dv/dt = 0, yields: 

K{to) fio\f awgsincto \ .„„s 

"0 - 777 '■ — ^' ^0 - ~ ■ n ^ u - ^0 (37) 

L{to,xo,xo,yo) \yo J \sm{l + c)to J 



where 
and 



K = bc^/^ sin(l + c)to 



L = ( 2xobc'/^ cos(l + c)to + % + ^^^^] + ahc^'^ sin etc 



The stationary point (Mo,fo) is a saddle (hereafter called 'X-point'), with one positive 
and one negative real eigenvalues. The reality of eigenvalues follows immediately by 
noticing that the variational matrix of (I25p is symmetric by virtue of the fact that the 
equations of motion are given by the grad Vu,vS' with S' = S — xqu — yov, with S{u,v, t) 
equal to the phase of the wavefunction ip = Re^^ . Thus, the off-diagonal elements of 
the variational matrix are equal, namely: 

d du d^S' d^S' d dv 



dv dt dvdu dudv du dt 

i.e., the variational matrix is symmetric and its eigenvalues are real. Furthermore, 
setting Oil = d{du/dt)/du, ai2 = d{du/dt)/dv, 021 = d{dv / dt) / du, 022 = d{dv/dt)/dv, 
the characteristic equation is given by 



— (an + a22)A + (011022 — 0210-12) = 
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To the lowest approximation we find: 




+ ..., 



Thus the roots of the characteristic equation satisfy 




/ dG OGV 
\ dv du j 



(38) 



and if we replace u, v in fl55]) by the root Uq, fo, with mq, Vq sufficiently small, the product 
of the eigenvalues is negative. This means that one eigenvalue is positive and the other 
negative (except for degenerate cases in which one eigenvalue is zero). 

A further conclusion stems by noticing that, if the X-point has a distance do = 



+ Vq from the nodal point, implying that both uq and Vq are of order 0{do), then 
all the entries Oij of the variational matrix are of order ©(I/c/q). It follows that both 
eigenvalues satisfy: 



This conclusion is important because it implies that while, as we will see in the next 
subsection, chaos is introduced mainly at the approch of the orbits near an X-point, the 
contribution of the latter to the positive value of the Lyapunov characteristic exponent 
of an orbit is determined by the measure of the X-point's positive eigenvalue A, which, 
on its turn, is large when the X-point is close to the nodal point, i.e., when d^ in Eq. fl39l) 
is small. This means that the nodal points influence chaos rather indirectly, that is, 
the chaotic behavior is actually due to the X-points, but the effectiveness of the latter 
depend on their closeness to the nodal points. Notice that the X-point can approach 
arbitrarily close to the nodal point, since the two points collide whenever K = 0, i.e., 
sin(l + c)to = 0. This happens whenever the nodal point reaches infinity from either 
side of the y-axis. In general, the distance do is small when \yo\ is large. 

Figure 6 shows how do the spirals emanating from the nodal point connect to the 
invariant manifolds emanating from the X-point in the adiabatic approximation. This 
figure is a numerical calculation of all the integral curves emanating from the X-point, 
when a = b = 1, c = and to = 10. The X-point is found numerically up to 

twelve digits by looking for roots of Eqs.f l25l) close to the nodal point by a Newton- 
Raphson routine. The numerical solution {uo,Vo) is then inserted in the expressions 
for the matrix elements aij,i = 1,2, j = 1,2, yielding the eigenvalues and eigenvectors 
of the variational matrix at {uo,Vo)- Then, we give initial conditions for u,v on both 
semi-lines (with respect to (mq, fo)) determined by the directions of the two eigenvectors, 
at a distance 10~^ from the X-point. Finally, we integrate numerically the differential 
equation: 





(39) 




(40) 



Nodal points and the transition from ordered to chaotic Bohmian trajectories 18 



0.00 



-0.50- 



1 .00+-r 

-0.50 




0.00 



0.50 



1 .0e 



u 



Figure 6. The nodal point - X point complex on the (u, v) plane, in the adiabatic 
approximation when = 10 and a = 6 = 1, c = \/2/2. 



found by dividing the two equations of (!25l) . for each of the four different above sets of 
initial conditions. This yields numerically the two branches of the stable and unstable 
manifolds of the X-point. Clearly, since all the spirals terminating at the nodal point 
are described in the same sense, only one of these four branches can be connected 
to a spiral terminating at the nodal point. This branch can always be identified by 
comparing the senses of description of the manifolds and of the spiral. In particular, 
one of the asymptotic spirals of the nodal point is joined to one branch of the unstable 
manifold emanating from the X-point, if the nodal point is an attractor, or to the stable 
manifold, if the nodal point is a repellor. The set of all the integral curves of the flow 
in the neighborhood of the nodal point and X-point is hereafter called the 'nodal point 
- X-point complex'. 

Figure 7 shows a comparison of the spatial distribution of the nodal points (Fig. 7a) 
and of the respective X-points (Fig. 7b), on the plane (x, y) when a = b = 1, c = \/2/2 
and to is in the interval < to ^ 1000. The nodal lines and the X-point lines form 
similar patterns. In particular, similarly to the nodal points (section 2), the X-points 
avoid a central region of the plane {x, y) in which the majority of orbits turn to be 
regular (subsection 4.3). Fig. 7c shows the modulus of the positive eigenvalue of the 
X-point as a function of the distance d^ of the X-point from the nodal point. We 
find numerically a power-law scaling A oc l/# with p ~ 1.5, i.e., less steep than the 
theoretical scaling given by Eq. (l39ll . From this figure, as well as from Figure 6, in which 
the distance of the X-point from the nodal point is about d^ = 0.9, we deduce that the 
results obtained by the previous perturbative analysis are essentially valid not only at 




Figure 7. (a) Nodal lines, (b) X-point lines, and (c) the positive eigenvalue of the 
X-point versus the distance do of the X-point from the nodal point, when a — b — 1, 
c = V2/2. 



very small distances from the nodal point but also at relatively large distances (of order 
10^^). At any rate, we always find that the stationary points of the vector field (!25|l . as 
determined numerically, induce a similar phase-space structure as in Figure 6, i.e., this 
structure is general in the model considered. 

4-2. The exponential sensitivity of the orbits 

In order to understand how does the approach of an orbit to the nodal point - X-point 
complex introduce exponential sensitivity of the orbits to the initial conditions, we 
consider in detail the successive encounters of two nearby orbits, with initial separation 
10~^, with this complex, which take place at snapshots at which the orbits pass close 
to the complex. To this end, we consider the orbit of Figure 8a (initial conditions 
xi(0) = 1/1(0) = —1.1 and a = b = 1, c = V2/2) which has a number of consecutive 
encounters with the nodal point - X point complex. This orbit is chaotic, as seen from 
the calculation of the 'finite time Lyapunov characteristic number' 

' t 1^(0)1 ^ ^ 

where ^{t) = {dx(t), dy{t)) is the deviation vector associated with the orbit (xi(t), ?/i(t)), 
which is found by integrating the variational equations of motion together with 
the orginal equations of motion. The limit lim^^oo x(^) yields the usual Lyapunov 
characteristic number. Numerically we find (Fig. 8b) that this limit is close to LCN ~ 
3 X 10~^. We then consider in detail the growth of deviations from this orbit by 
calculating also a nearby orbit {x2{t),y2{t)) with initial conditions ^2(0) = a;i(0) + 10~^, 
7/2(0) = 1/1(0). 

Figure 9a shows the growth in time of the distance AS{t) = ( (xi(t) — X2{t))'^+ 



\ 1/2 

{yi{t) — 7/2 (^))^ I between the two nearby orbits. Clearly, the distance grows in 
general with time, but the growth takes place by rather abrupt steps. That is, while the 
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Figure 8. (a) Chaotic orbit with initial conditions x(0) = y{0) = —1.1, and a = b = 1, 

c = •\/2/2. The orbit intersects the nodal lines of the left semiplane of {x,y). (b) The 
'finite time' Lyapunov characteristic number x{t) for the same orbit. The limiting 
value is close to LCN ~ 0.03. 

distance has in general large fluctuations, there are particular times when the distance 
AS suddenly grows by jumps of about one order of magnitude (or more). Thus, the 
initial distance AS — 10~^ becomes of order 10~^ at a time t ~ 25, then of order 10"^ 
at i ~ 90, and finally of order 10"^ (reaching even unity) at about t ~ 170. After this 
time the distance A^" can no longer be considered as small, that is, the orbits are no 
longer nearby. 

Figure 9b is a close up to the third of the above described jumps, focusing on the 
time interval 160 < t < 180. From this figure it is clear that there are three encounters 
of the orbit with the nodal point - X point complex taking place in the considered time 
interval. In encounter (1), the minimum distance of the orbit (1) (= {xi{t),yi{t)) from 
the nodal point is emin — 0.4, and the minimum distance from the X point is even smaller 
(d^in — 0.1). On the contrary, in the next encounter (II), the minimum distance from 
the X-point is rather large {dmin — 0.9) while the minimum distance from the nodal 
point is about the same as in case (I). Finally, in case (III) the minimum distance from 
the X-point is small {dmin — 0.2) while the minimum distance from the nodal point is 
now large (e„jj„ ~ 1). Clearly, the growth of the distance AS mostly takes place during 
the encounters (1) and (III) in which the orbits pass closer to the X-point than to the 
nodal point. On the other hand, in the case of encounter (II), the growth is smaller 
while the orbit approaches closer the nodal point than the X-point. We conclude that 
large variations of AS are in general associated with approaches of the orbits to the 
X-point rather than to the nodal point. 

Figure 10 shows in detail how does the separation of the orbits take place during 
the encounter event (III). The two nearby orbits are shown as dashed curves, in the time 
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Figure 9. (a) The separation AS of two nearby orbits as a function of the time t 
{a = b = 1, c = and initial conditions a;i(0) = —1.1, a;2(0) = a;i(0) + 10"'*, 

yi(0) = t/2(0) = -1.1). (b) A detail of (a) in the time interval 160 < t < 180. The line 
with thick dots gives A5. The solid and dashed lines show the distances d{t) and e{t) 
of the orbit {xi{t),yi{t)) from the instantaneous locations of the X-point and of the 
nodal point respectively, when the latter are smaller than 1. 

interval 175 < t < 176.5, in the moving frame of reference centered at the nodal point. 
The different frames correspond to different time snapshots, and the particular positions 
of the orbital points (1) = {ui, Vi) = {xi - xq, yi - yo), (2) = (^2, V2) = (^2 - xq, 2/2 - yo) 
at the given snapshot are marked with thick dots. Finally, the instantaneous stable 
and unstable asymptotic manifolds emanating from the X-point (in the adiabatic 
approximation) are plotted for the time corresponding to each frame. 

We notice that the X-point changes position relative to the nodal point, which in 
these frames is always centered at {u,v) = (0,0). In fact, as already mentioned, the two 
points collide whenever sin(l + c)t = (and then yo = ±00). There are two consecutive 
collisions at t = 174.829 and t = 176.669. The time t = 175.2 (Fig. 10a) is close to 
the first of the above two collision times and, consequently, the X-point at t — 175.2 
is very close to the nodal point. Then, the X-point moves to the left up to about 
t — 175.8 (Figs.lOb,c,d), and then it returns to the right approaching again the nodal 
point (Figs.l0e,f). The time t — 176.3 is close to the second collision time {t — 176.669), 
thus the X-point in Fig.lOf comes again very close to the nodal point. 

Now, at the time t = 175.2 the orbits have a separation of about AS* = 0.15 
(Fig. 10a). At this time snapshot the orbits move in a nearly parallel way, and their 
distances from both the nodal point and X-point are rather large (of order unity). The 
orbits move downwards in about the same direction as indicated by the arrows of the 
invariant manifolds of the X-point (in the rest frame {x, y) this means that the nodal 
point approaches the orbital points from y — —00, see also Fig. 11a). Furthermore, as the 
X-point itself moves from right to left, both orbits approach to it (Fig. 10b, t — 175.5). 





Figure 10. The encounter event (III) of Figure 9b viewed in detail in the plane v). 
The orbits (l)and (2), in the time interval 175 < t < 176.5, are shown by dashed lines, 
while the thick dots indicate the positions of the orbital points on these lines at the 
times (a) t = 175.2, (b) t = 175.5, (c) t = 175.7. (d) t = 175.8, (e) t = 176, and (f) 
t = 176.3. The invariant manifolds of the instantaneous X-point - nodal point complex 
are plotted for the same times. The main deflection of the orbits, takes place within 
the time interval 175.7 < t < 175.8. 
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Figure 11. (a) The nodal point coordinate yo(i) in the time interval 175 < t < 176.5. 
(b) The value of < /a > (to) in the time interval 175 < to < 176.5. There are three 
moments in this interval at which < /s >= 0. These correspond to changes in the 
topological structure of the X-point - nodal point complex. 



The crucial phenomenon occurs near t — 175.7 (Fig. 10c). Around this time, the 
X point crosses a segment joining the orbital points (1) and (2). The two points have 
approached the moving X-point at a distance smaller than 0.2, but they are on opposite 
branches of the unstable manifold. Thus, point (1) moves downwards following one 
branch of the unstable manifold of the X point, while point (2) moves upwards following 
the other branch of the same manifold. This causes a abrupt growth of the distance of 
the two points by a factor ~ 3. 

An important change in the topological structure of the invariant manifolds, that 
influences the orbits, takes place between t — 175.7 (Fig. 10c) and t — 175.8 (Fig.lOd). 
Namely, a,t t — 175.7 (Fig. 10c) the spiral terminating at the nodal point is connected 
to one branch of the stable manifold of the X-point, while at t = 175.8 (Fig.lOd) it is 
connected to one branch of the unstable manifold of the X-point. This transition takes 
place via a Hopf bifurcation which is described in detail below. At any rate, at t = 175.8 
(Fig.lOd) the X-point has moved to the left, far from the nodal point, and orbit (1) is 
close to the stable manifold of the X-point. Thus, orbit (1) is deflected to the left, while 
orbit (2) continues slowly upwards. 

Finally, a little later {t — 176, Fig.lOe), the X-point returns close to the nodal point 
so that point (2) comes very close to the unstable manifold of the X-point. This causes 
a deflection of orbit (2) to the left, while orbit (1), although far from the X-point - 
nodal point complex, moves also to the left, following the general direction of motion 
indicated by the unstable manifold of the X-point. Finally, at t = 176.3 (Fig.lOf) both 
orbital points are far from the X-point - nodal point complex, but they are relatively 
close to the unstable manifold of the X-point in a downwards direction. From there on 
both orbits move in a nearly parallel way until the next encounter event which occurs 
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Figure 12. The instantaneous flow cliart of ttie X-point - nodal point complex at 
the times (a) to = 176.7, (b) t = 176.76, (c) t = 176.7750938, (d) t = 175.78. A Hopf 

bifurcation taking place near the value t = 176.72 (at which < >= in figure 11), 
leads to the formation of a limit circle (repellor), shown in (b), which disappears after 
reaching the separatrix of (c). 

much later. The overall growth of the distance of the two orbits by a factor 3 in a 
time At — 176.3 — 175.2 = 1.1 corresponds to an exponential growth rate In 3/1.1 ~ 1 
in this time interval. This is much larger than the average exponential growth rate 
(=LCN~ 0.03) for the same orbit within a much longer time interval. This fact justifies 
the statement that the growth of AS is by abrupt jumps, which take place during local 
(in space and time) encounters with the nodal point - X-point complex. 

The topological transition in the phase space structure taking place between 
t — 175.7 (Fig. 10c) and t = 175.8 (Fig.lOd) is due to a Hopf bifurcation taking place 
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between these two times as a result of the fact that the value of < /s > (Eg. (1551) ) 
changes sign (Fig. lib), while the sign of d(f)/dt remains constant (because there is no 
change in the sign of sin(l + c)t, i.e., no transition of yo to infinity, in the same interval, 
Fig. 11a). By virtue of Eq. (l36|) . the change of the sign of < /s > at a time between 
to = 175.7 and to = 175.8 implies that the nodal point turns from repellor to attractor. 
The precise time when this happens depends on higher order terms in the development 
of the equations of motion around the nodal point, but it is nevertheless close to the time 
when the term depending on < /s > becomes equal to zero. Numerically, we find the 
bifurcation to take place near to = 175.75. Before this time (e.g. at to = 175.7, Fig.l2a) 
the nodal point is a repellor, and a spiral emanating from it joins the stable manifold 
of the X-point. On the other hand, a little after this time (to = 175.76, Fig. 12b), the 
nodal point has become an attractor, while the stability character of the X-point does 
not change appreciably. Thus, between the nodal point and X-point there is now a limit 
circle which acts as a repellor, i.e., the orbits on both sides of the circle move away from 
it on spirals either terminating at the nodal point or moving towards the X-point. In 
the latter case a spiral either joins the stable manifold terminating at the X-point or 
continues downwards, away from the X-point, in the channel formed between the two 
branches of the unstable manifold of the X-point (Fig. 12b). As to increases the limit 
circle moves outwards approaching the invariant manifolds of the X-point. At a critical 
time to = 175.7750938 (Fig. 12c) the limit circle coincides with the invariant manifolds of 
the X-point. At still larger times (t = 175.78, Fig.l2d), the limit circle disappears and 
one branch of the unstable manifold of the X-point continues now as a spiral terminating 
at the nodal point. 

Figure 13 shows how does the effect of a close encounter of an orbit with the 
nodal point - X point complex shows up in the time evolution of the deviations ^(t) 
as given by solving the variational equations of motion together with the equations 
of motion for one orbit. The abrupt jumps in the length of the deviation vector 
^ = {dx,dy) = {du,dv) (Fig. 13a) are associated with passages of the orbit close to 
the nodal point - X-point complex. In order to obtain quantitative estimates of the 
exponential growth of deviations at successive encounter events, we proceed as follows: 
The time span of the total run of an orbit is split in short windows of width At = 0.1 (the 
timestep of the numerical integration is variable and much shorter than this window, 
i.e., dt < 10^^). In each time window, an average 'stretching number' (Voglis and 
Contopoulos 1994) is calculated according to: 

1 , ,f(ti + At), , , 

At ' ' ^ ^ 

where tj = (z — 1) x At is the initial time of the i-th window, and ^(t) is the length of the 
deviation vector ^ (t) at the time t. This quantity characterizes the local growth rate of 
deviations, while the average value of all the stretching numbers yields the 'finite time 
Lyapunov characteristic number' 

1 " i ^fe+Ai) 1 m, 
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(The limit limt_>oo x(^) defines the Lyapunov characteristic number of the orbit). In the 
same time windows we also store the values of the minimum distance of the orbit to the 
nodal point emin, and to the X-point dmin, as well as the minimum distance between 
the X-point and nodal point do^rnin- These three minimum values are not occurring at 
precisely the same times within a given time window, however the width of the window 
At — 0.1 itself is small enough so that the occurrences of the minimum values can be 
considered as nearly simultaneous. 

Figure 13b then shows the main result. The distances that an orbit reaches from 
either the nodal point or the X-point, for all the 10000 time windows in a total time 
span < t < 1000, are grouped in bins of width S — 2.5x 10~^. The abscissa in Fig. 13b 
gives the median value of each bin, which is the same for either the distance d from the 
X-point or e from the nodal point. A weighted average value a of the stretching number 
is then calculated in each bin, by summing the values of the stretching numbers which 
appear during all the passages of the orbit at distances d ± 5/2 from the X-point, or 
e±5/2 from the nodal point, and dividing by the total number of values of the stretching 
number in the sample. Clearly, when the stretching numbers are grouped with respect 
to the various distances d reached from the X-point, the average stretching number 
a{d) is positive for all distances d < 0.25, while for d > 0.25, d{d) fluctuates between 
positive and negative values, showing nevertheless a preference for positive values. Such 
a preference reflects the hyperbolic dynamics induced on the orbits by their approaches 
to the X-point. That is, the general solution of the variational equations close to the 
X-point contains terms growing exponentially and other terms decaying exponentially. 
However, the growing terms prevail as the time t increases. Thus, the overall average 
stretching number < a > after many encounters of an orbit with the X-point turns to 
be positive. 

On the other hand, when the stretching numbers are grouped with respect to the 
distance e of an orbit from the nodal point, the average stretching number a(e) (dashed 
curve in Fig. 13b) shows no clear preference towards positive or negative values when e is 
small (e < 0.25), while a preference towards positive values of a appears when e > 0.25. 
This agrees with Fig. 13a, or Fig.9b, which show that the growth of deviations occurs 
mainly during encounters (1) and (111), during which the minimum distance of the orbit 
from the X-point is smaller than the distance from the nodal point. 

5. Discussion 

The role of the nodal points of the wavefunction (sometimes called 'quantum vortices') 
in determining the main features of the 'hydrodynamical' probability flow has been 
studied in a number of different quantum systems (see Wyatt (2005) and references 
therein). Here we refer only to studies related to our own, i.e., to the appearance of 
chaos due to nodal points (quantum vortices). 

Berry (2005) studied the flow hues of a general time-independent complex scalar 
field ip{x,y) and found that they typically spiral in or out of a 'phase vortex'. He also 
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Figure 13. (a) the growth of the length of the deviation vector £^{t) for the same 
orbit {xi{t),yi{t)) as in Figure 9b. The lower solid and dashed lines are the same 
as in Fig. 9b. (b) The mean stretching number a as a function of the distance from 
the X-point a{d), or from the nodal point a(e), for the same orbit and time up to 
t = 1000. For each value of d or e, the value of a is calculated by taking the sum of all 
the stretching numbers (Eq. ([^ with At = 0.1) occuring at passages of the orbit at a 
distance dztS/2 from the X-point and e ± 6/2 from the nodal point respectively, with 
5 = 2.5 X 10-2. 



found the stationary points of this flow which can be either elhptic or hyperbohc. Berry 
exphcitly excludes a spiral flow near the vortices of the Schrodinger field because the 
quantum mechanical current satisfies V ■ j = 0. This is precisely what happens in our 
case if one considers the flow lines in the rest frame {x, y) rather than a frame moving 
together with the nodal points. In the rest frame, the flow integral curves (found by 
dividing by parts Eqs.(jl])) are given by the differential equation: 

dy bc^^'^x{ax sin ct + sin(l + c)t) 



dx a sin t + hc^/'^y sin(l + c)t 
Keeping the time t frozen in the r.h.s., Eq. fH^ yields the integral curves 



(44) 



2 1 

a sin ty + 

6c^/2 sin(i + cy— - ahc^'^ sin ct— - hc^'^ sin(l + c)t— = C(45) 

The critical points of (14 5 p are given by the solutions of VC = 0. There are two solutions: 
, , , sin(l + c)t asint , , , 

^"^'^^^ = ^'^;i^'- &cV^sin(l + c)t ^ " 

The first critical point (xi,?/i) coincides with the nodal point. The eigenvalues of the 
Hessian matrix of C at (a;o,?/o) are given by Ai^2 = ±i&c^/^ sin(l + c)t, thus they are 
imaginary at any time t, implying that the integral curves of the velocity field in the 
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neighborhood of the nodal point are approximately ellipses centered at (xo,?/o) at any- 
time t except t(. = 2k7c/{l + c) or tc = {2k + 1)tt/{1 + c) with k integer. The difference 
with respect to the approximation of Eq.( l36l) is that in the moving frame of reference 
the instantaneous flow lines form spirals if < /s 0, i.e., if io 7^ or ?/o 7^ 0, that is 
the spirals appear (in the moving frame) only because the velocity of the nodal point 
is non-zero. Furthermore, in our analysis, the eigenvalues of the X-point in the nodal 
point - X-point complex scale as an inverse power of the distance of the X-point from 
the nodal point. On the contrary, the eigenvalues of the second critical point {x2,y2) of 
Eqs.( H6|l are given by Ai^2 = ±&c^/^sin(l + c)t, i.e. they are bounded by quantities of 
order a or b. Thus, despite the fact that (x2, ^2) represents also a saddle point in the rest 
frame of motion in the adiabatic approximation, its influence to the dynamics is not so 
important when a and b are of order unity. In fact, this point is always attached to the 
y-axis, so that it can only influence the deviation vectors at times when the orbits come 
close to this axis. This is at variance with the numerical results showing that positive 
stretching numbers are introduced when the orbits approach the nodal point - X-point 
complex at arbitrary locations within the plane {x,y). 

Wisniacki and Pujals (2005) studied the Bohmian orbits in an example similar 
to ours, namely the superposition of three stationary states of the double harmonic 
oscilator, when, however, the latter is isotropic (c = 1) and the ratios of the probability 
amplitudes of the states |\['io > and |\E'ii > with that of |\E'oo > are complex. In the 
case 0=1 there is only one fundamental frequency of the time- dependent trigonometric 
terms of the equations of motions. This allows one to obtain stroboscopic plots of 
the orbits, i.e., Poincare surfaces of section. Wisniacki et al. pointed out that it is 
the motion of the nodal point which generates chaos. However, their mechanism of 
introduction of chaos is different from ours. Namely, in the case of Wisniacki et al. the 
motion of the nodal point generates a saddle point on the surface of section, but the 
surface of section is area preserving, and the transverse intersections between the stable 
and unstable invariant manifolds of the saddle point generate homoclinic chaos. The 
most important difference is that their mechanism can be applied only in resonant cases 
while our mechanism applies to general non-resonant cases. 

A similar example was studied by Makowski et al. (2000). Both results correspond 
to the case of a complex ratio a/b and c = 1. On the other hand, if c is rational but a/b 
is real, all the orbits are periodic and neutrally stable, thus there is no chaos at all (some 
examples of seemingly chaotic orbits in a similar model given by Konkel and Makowski 
(1998), for a rational value of c, are just due to numerical errors caused by the stiffness 
of the equations of motion close to the X-point. In order to obtain the correct orbits, 
which are periodic and not chaotic, we had to use a program in Mathematica with an 
accuracy of 50 digits!). 

Falsaperla and Fonte (2003) studied the orbits near nodal lines in a 3D model. 
In that case the orbits describe helical motions around the nodal line (called by these 
authors 'spirals') while the projections of the motion on a z =constant surface are 
ellipses. These authors point out that the nodal lines "regularize the motion" in their 
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neighborhood, and only orbits not following a definite nodal line are "intermittent 
chaotic". They furthermore find too that the period of rotation along the helix is 
of order O(e^), where e in this case is the distance from the nodal line. While a detailed 
comparison of 2D and 3D models is necessary, our analysis above shows essentially why 
the growth of the deviation vectors is intermittent, i.e., it takes place by abrupt steps 
whenever an orbit approaches the nodal point - X-point complex. Furthermore, we also 
explain why the appearance of chaos is not strictly correlated with very close approaches 
to the nodal point (in such cases the orbits simply spiral around the nodal point) , but it 
occurs mainly when the orbit approaches closely the X-point associated with the nodal 
point. This happens even if the latter has a distance from the nodal point which is much 
larger than the distance of the orbit from the X-point. 

Wu and Sprung (1999) gave plots (their Fig. 4) of the probability flow, which 
resemble our Figure 7. However, this resemblence is only due to a 'phase mixing' 
phenomenon taking place in the rest frame near the nodal point. Namely, because 
the rotation frequency depends on the distance from the nodal point as a fluid 

element of some thickness approaching the nodal point forms a number of windings 
around the nodal point due to the differential rotation of the orbits included in its area. 
These windings give the impression of forming a spiral pattern, which is, however, only 
apparent, namely the windings arc limited by an inner circle (or ellipse) due to orbits 
of the fluid element closest to the nodal point (this is different from our limit circle of 
Fig. 12). On the other hand, in earlier works the same authors found spiral motions 
around quantum vortices in the rest frame of motion, when the quantum system is 
subject to a vector potential in addition to a scalar potential (Wu and Sprung 1994). The 
plots shown in that work are similar to ours, although the similarity in the topological 
structure of the phase space is due to the presence of an extra term in the equations of 
motion due to the vector potential, rather than to the motion of the nodal point. 

Finally, Frisk (1997) and Wisniacki et al. (2006) noticed that the transition to 
chaos is enhanced when there are many co-existing vortices inffuencing the orbits. In 
that case it is necessary to consider the connections of the invariant manifolds emanating 
from the stationary points of one vortex with those of other vortices and see how this 
can increase chaos and transport phenomena in conflguration space. This problem is 
proposed for future study. 

6. Conclusions 

In summary, the main conclusions of our study are the following: 

1) In a simple quantum system consisting of the superposition of three eigenstatcs 
of the double harmonic oscillator potential, in which there is one nodal point of the 
wavefunction travelling in the plane {x,y), we proved the existence of domains of this 
plane which are free of nodal points. Bohmian orbits in these domains, as well as orbits 
slightly overlapping with the nodal lines, are regular. 

2) In the central domain devoid of nodal lines the equations of motion admit 
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expansions in powers of the probability amplitudes of the eigenstates, yielding the 
trajectories as double Fourier series in the fundamental frequencies of the system. 
These series represent theoretical orbits which are, by definition, regular. We show 
the agreement of the theoretical and numerical orbits in this domain for sufficiently 
small values of the probability amplitudes and sufficiently high order of the expansions, 
and use this to explain the morphological characteristics of the regular orbits. 

3) Close to the nodal points, we use different expansions, in powers of the distance 
from the nodal point, in order to unravel the dynamics. The angular frequency of motion 
in a e— neighborhood of the nodal point is of order 0(l/e^). This justifies the use of 
the adiabatic approximation. In the fixed frame {x, y) the flow lines close to the nodal 
point are ellipses. However, in a moving frame attached to the nodal point, the flow 
lines are spirals terminating at the nodal point. The temporary sense of description 
of the spirals by the orbits is unique, i.e., at a given time the nodal point is either an 
attractor or a repellor. Furthermore, at a flnite distance from the nodal point there is a 
saddle stationary point of the flow with one asymptotic manifold joining the spiral and 
the other three extending to inflnity. 

4) The eigenvalues of the X-point scale as I/^q where do is the distance of the 
X-point from the nodal point and p ~ 2. Furthermore, the distance do can be arbi- 
trarily small, i.e., there are collisions of the X- and nodal points. We show that the 
orbits approaching close to the X-point exhibit exponential growth of their deviations, 
i.e., they arc chaotic. In all numerical examples we find that chaos is associated with 
the approach of the orbits to the X-point, which, however, is only guaranteed when we 
have a moving nodal point. Furthermore the chaotic infiuence of the X-point on the 
orbits is strongest when the X-point is closest to the nodal point. Thus, the nodal point 
indirectly influences the transition of Bohm's trajectories from order to chaos via the 
above mechanism. 
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7. Appendix: Calculation of < > 



The precise form of Eq.f l30p reads: 

dR —AR^ — xqG cos (j) — yoG sin i 



d(j) B — CR — cos (p + Xq^ sin (p 

where 

A = abc^^'^ cos"^ (f) sin (f) sin at, B = bc^^"^ sin{l + c)t, C = abc^^"^ cos'^ (j) sin at 
and G = g2R^ + g^R^ + g^R^ with 



(47) 



g2 = — ir~ ~ 26c^/^ sin cos cos(l + c)t + 6^cXq sin^ > 



.2 



(73 = — — — cos(l + c)t cos^ sin + 2b^cxQ cos sin^ 

(74 = cos^ sin^ . 
Expanding Eq. (l47l) in ascending powers of R we find: 

^ = ^ - AR^ - xog2 cos (pR^ - x^g^ cos 0i?^ - yQg2 sin (pR^ - yogs sin 0i?^ j 
X -^h +-|i? + ?7o§cos0i?-i;o§sin0i?j + ... (48) 



dR I r^^ dR , . 9 . ^ 

-rf0=</2>i?2+</3>i?' + ... (49) 



Averaging Eg. (1481) we find: 

-/ 

d(p 2tc Jo 
where ^ 

/a = ^(-^ - io5'2cos0 - ?/ofi'2sin0) 

/a = -^(-^-a;ofl'2cos0-?7ofl'2sin0) x (-^ + cos - sin , 

+ ;^(-a;ofl'3Cos0 - ?7ofl'3sin0) . 

Collecting terms of the form 

1 /•2^ 
27r io 

we find that < /a >= 0, while 

< /3 > = ,652, ^ (^-bc^^ ^inct(-| + 6^ccos(l + c)t^o) 

+ 4xoyo(&'c2x^ - A) - 4&ci^' cos(l + c)t(i;2 _ yo')(A + b^cxD] 

Xq Xq ) 

I ^C^^' / j/oCOs(l + c)to , l/2^.\ ..^N 

+ 46cV2sin(l + c)tol ° V ^ ^ 

which after some simplification yields Eq. (l35|l . 
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